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This paper reports on our efltort in modeling realistic astropiiysical neutron star binaries in gen- 
eral relativity. We analyze under what conditions the conformally flat quasiequilibrium (CFQE) 
approach can generate "astrophysically relevant" initial data, by developing an analysis that de- 
termines the violation of the CFQE approximation in the evolution of the binary described by the 
full Einstein theory. We show that the CFQE assumptions significantly violate the Einstein field 
equations for corotating neutron stars at orbital separations nearly double that of the innermost 
stable circular orbit (ISCO) separation, thus calling into question the astrophysical relevance of the 
ISCO determined in the CFQE approach. 

With the need to start numerical simulations at large orbital separation in mind, we push for 
stable and long term integrations of the full Einstein equations for the binary neutron star system. 
We demonstrate the stability of our numerical treatment and analyze the stringent requirements on 
resolution and size of the computational domain for an accurate simulation of the system. 
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I. INTRODUCTION 

The analysis of general relativistic binary neutron star 
processes is an important, yet challenging, endeavor. The 
importance of understanding these processes is rooted in 
observational astronomy, both in gravitational wave as- 
tronomy and high-energy electromagnetic wave astron- 
omy. Neutron star binaries could be the central en- 
gines for some classes of gamma-ray bursts, and they are 
definitely strong candidates as sources of gravitational 
radiation detectable by the up and coming generation 
of gravitational wave detectors such as LIGO, VIRGO, 
TAMA, GEO, and LISA. The challenge of understanding 
binary neutron star processes is rooted in the complex- 
ity of both the nonlinear Einstein field equations and the 
physical properties of the super-nuclear density matter 
which make up the neutron stars. 

Due to the complexity of the binary neutron star sys- 
tem, various levels of approximation have been employed 
(with various levels of success) as aids in understanding 
the details of the different stages of the inspiral of compa- 
rable mass binary neutron stars, from the quasistationary 
inspiral stage through to the plunge and merger of the 
binary stars to the ring-down of the final merged object. 
These range from the approximation of the structure of 



the neutron stars themselves (e.g., from point particle 
to finite sized perfect fluid models and the equations of 
state with different physical approximations) to the ap- 
proximation of general relativistic effects (e.g., from New- 
tonian gravity, to post-Newtonian, to full general relativ- 
ity). Based on insights from point particle mechanics in 
general relativity and orbits of finite size bodies in Newto- 
nian gravity, one expects that the early part of the inspi- 
ral process will be quasistationary, with the secular (i.e. 
on timescales larger than the orbital timescale) shrinkage 
of the orbit driven by gravitational radiation. When the 
orbit shrinks to a small enough radius (but before touch- 
ing), there may or may not exist an innermost stable 
orbit (ISO) beyond which dynamical processes drive the 
quasistationary inspiral into a plunge phase. Whether 
there is, in fact, a "phase change" from quasistationary 
inspiral to a plunge phase in the fully relativistic the- 
ory, and whether or not this happens before any other 
(hydrodynamical) dynamical instabilities, is an unsolved 
problem in full general relativity (see, e.g., P| for answers 
in the Newtonian case). 

A recent approach for the investigation of the later part 
of the neutron star inspiral which has been studied in de- 
tail 2, 3, 4, 5, 6, 7, 8j 9,_lQa JJJ has drawn much attention. 
This treatment, which we refer to as the conformally flat 
quasiequilibrium (CFQE) approach, is a procedure for 
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constructing general relativistic configurations that cor- 
respond to compact, equal mass binary neutron stars in 
a quasiequilibrium, circular orbit. These individual con- 
figurations, which we refer to as "CFQE configurations" , 
are by themselves solutions to the constraints of general 
relativity in 3+1 form, i.e., the Hamiltonian and momen- 
tum constraints. For a given equation of state, each equal 
mass binary CFQE configuration can be characterized by 
two parameters: the separation of the two neutron stars 
and the baryonic (rest) mass of each of the neutron stars 
(we only consider equal mass, corotating binary configu- 
rations here) . Taking the CFQE approximation one step 
further, one can construct an entire 4-dimensional space- 
time by "gluing" these configurations together as a time 
sequence of different CFQE configurations. This con- 
struction, which we will refer to as a CFQE-sequence, 
produces a spacetime that solves 5 of the Einstein field 
equations (4 constraint equations and trace of the ex- 
trinsic curvature equation) 0, . Since the total bary- 
onic mass will be constant during the inspiral, the CFQE 
approach takes constant baryonic CFQE-sequences as a 
representation of the evolutionary sequence of the binary 
neutron star inspiral process; secular evolution of the or- 
bit due to gravitational radiation brings one equilibrium 
configuration into another with the same rest mass, form- 
ing the CFQE-sequence 0, H ES El . 

In this CFQE-sequence approach, one finds that for 
a corotational neutron star binary system (which is the 
focus of this paper), the ADM mass of each of the indi- 
vidual CFQE configurations decreases along a constant 
rest mass CFQE-sequence as the separation between the 
binary stars decreases until a minimum is attained, af- 
ter which the ADM mass increases as the separation de- 
creases further. Using turning point theorems for exact 
equilibrium configurations in general relativity [3. ll2l[T^ . 
this minimum point signals a secular instability in the 
evolutionary sequence, and is commonly referred to as 
the innermost stable circular orbit (ISCO) configuration. 
In the Newtonian case, it has been shown that this sec- 
ular instability point is encountered in the evolutionary 
process before any dynamical instability is reached |l4j . 

While it is reasonable to assume that the CFQE ap- 
proximation holds reasonably well for highly separated 
neutron stars, it is not clear at which point along the 
evolutionary sequence the CFQE approximation breaks 
down. It is certainly not clear whether or not the CFQE- 
sequence approximation is good for all neutron star sep- 
arations larger than the ISCO separation, due to the 
fact that one main assumption in the turning-point theo- 
rems ^.jAi Jij used to interpret these CFQE sequences is 
the assumption of exact equilibrium. While CFQE con- 
figurations are very close to equilibrium for large neutron 
star separations, they become less so as the separation de- 
creases. Unfortunately, it is exactly the small separation 
regime, i.e. near the ISCO configuration, in which these 
theorems are being applied. Also, the CFQE-sequence 
approximation to full general relativity does not provide 
an estimate of the error of its solutions, unlike, e.g., the 



post-Newtonian approach where one can compute the 
next order post-Newtonian terms to estimate the error 
of any post-Newtonian approximation. 

The fact that each CFQE configuration satisfies the 
constraints of general relativity suggests the use of these 
configurations as initial data to full general relativistic 
calculations. The setting of initial data is obviously an 
important issue in numerical general relativistic astro- 
physical simulations. While all initial data configurations 
satisfying the constraints of general relativity (i.e., the 
Hamiltonian and momentum constraints) are in principle 
legitimate initial data sets, they may not be acceptable 
for the study of coalescence of astrophysical neutron star 
binaries. In order for the results of the numerical evo- 
lutions to be relevant to observations, e.g., the gravita- 
tional waves emitted in actual neutron star coalescences, 
we have to make sure that the initial data actually cor- 
responds to a configuration in an astrophysically real- 
istic inspiral. It is not straightforward how one can go 
about evaluating the usefulness of the CFQE approach in 
approximating astrophysically relevant phenomena. We 
note that while the CFQE approach leads to solutions of 
the constraint equations, the CFQE-sequence approach 
is not consistent with the full set of Einstein equations. 

This paper is divided into four main sections. In sec- 
tional we describe our fully consistent general relativis- 
tic hydrodynamics code used in this paper. We describe 
in detail the 3-1-1 formulation of the Einstein field equa- 
tions that we couple to the relativistic hydrodynamics 
equations, along with details of the gauge conditions and 
discretization techniques used. 

In section lini we describe the CFQE-sequence approx- 
imation, whose individual configurations will be used as 
initial data for our fully consistent general relativistic 
calculations. We demonstrate that the CFQE-sequence 
definition of the ISCO for corotating neutron star bina- 
ries may not be a relevant concept by showing that if one 
takes into account the spin energy of the neutron stars 
when constructing the effective binding energy in the 
CFQE-sequence approximation, then there is no longer 
a minimum in the effective binding energy. 

In section IIVI we analyze the key assumptions of the 
CFQE-sequence approximation by comparing them with 
fully general relativistic simulations using CFQE config- 
urations as initial data. We focus on the question of 
how well the CFQE-sequence approximation for coro- 
tating neutron stars approximates full general relativity. 
This is done by comparing the assumptions in building 
the CFQE-sequence to fully consistent general relativis- 
tic calculations using CFQE configurations as initial data 
to our numerical evolution code. Specifically, we analyze 
the conformal fiatness assumption and the assumption 
of the existence of a Killing vector field. We devise a 
number of invariant measures of these assumptions, and 
monitor them in our fully consistent, general relativis- 
tic simulations. We find that, as expected, the accuracy 
of the CFQE-sequence approximation increases for in- 
creasing neutron star separation. We present a general 
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algorithm for evaluating whether or not any CFQE con- 
figuration can be thought of as astrophysically realistic 
initial data by analyzing how much the CFQE-sequence 
approximation violates the Einstein field equations. We 
demonstrate this method by showing that, for any given 
tolerance, one can find a CFQE configuration whose sub- 
sequent evolution in full general relativity will not violate 
the Einstein field equations within a small fraction of an 
orbit. 

In section [V] we analyze the long timescale (i.e. mul- 
tiple orbits of the binary system) numerical evolutions 
of CFQE initial data configurations using our general 
relativistic hydrodynamics code. In particular, we use 
a CFQE configuration which has a 15% larger proper 
separation than that of the ISCO CFQE configuration. 
We analyze the orbital decay rate on multiple orbit 
timescales. Using care in estimating the numerical trun- 
cation error as well as the errors introduced by placing 
the computational boundaries at a finite distance from 
the neutron star binaries, we find that while the com- 
putational resources at our disposal are sufficient to get 
a reasonable handle on the truncation error, the errors 
introduced by the placement of the boundary of the com- 
putational domain can have a significant impact on the 
dynamics of the neutron star inspiral. We demonstrate 
that this is true even if the linear dimensions of the com- 
putational domain are as large as one half the size of 
the gravitational waves being emitted, which is a fairly 
large computational domain by today's numerical rela- 
tivity standards. We conjecture that presently available 
computational resources will not allow a unigrid finite dif- 
ference code to decrease the discretization parameter Ax 
sufficiently and simultaneously increase the distance from 
the binary system to the computational domain bound- 
ary sufficiently in order to guarantee that the induced 
numerical errors will not significantly affect the details of 
the inspiral process. Mesh refinement techniques and/or 
better outer boundary conditions will be needed in order 
to accurately simulate the physics on the length scales 
of the compact object as well as the length scales of the 
gravitational radiation. 

II. FORMULATION AND DISCRETIZATION 
OF THE EINSTEIN EQUATIONS COUPLED TO 
A PERFECT FLUID 

Our code numerically solves the Einstein field equa- 
tions coupled to a relativistic perfect fluid. The gravi- 
tational degrees of freedom are geometrically encoded in 
the 4-metric, g^i,, which are governed by the Einstein 
field equations 

G^i, — SttT^^, (1) 

where G^i, is the Einstein tensor and T^j, is the stress- 
energy tensor of the perfect fluid, given as 

T^v = phu^u^ + Pgf_cu- (2) 



Here, we have set the gravitational constant G and the 
speed of light c to be identically 1. The 4- velocity of the 
perfect fluid is denoted as u^, and p, P, and h are the 
baryonic mass density, pressure, and specific enthalpy, re- 
spectively, of the fluid. The equations of motion govern- 
ing the perfect fluid are the conservation of stress-energy 
and baryonic mass 

y^j^t^i^ = 

V^.ipu'^) - 0. (3) 

These represent five equations governing the five degrees 
of freedom of the perfect fluid (the mass density, energy 
density, and velocity). The entire system of equations is 
closed by choosing an equation of state for the pressure 
P as a function of the baryonic mass density and internal 
energy density of the fluid. 

A. The Einstein equations in 3+1 form 

In order to numerically solve the Einstein fleld equa- 
tions, Eq. ^ we must cast the equations as an initial 
value problem. In order to facilitate this, we introduce 
a foliation of the spacetime into spacelike hypersurfaces 
where the coordinate t labels each spacial hypersurface. 
We furthermore introduce Cartesian coordinates on 
each spacelike hypersurface. The line element can now 
be written as 

ds^ = (/32 - a'^)dt^ + 2/3, dt dx^ + -fij dx' dx^ (4) 

where the shift vector /3* is a 3-vector defined on each 
spacelike hypersurface, a is the lapse function, and jij is 
the 3-metric. We denote 7'^ as the inverse of the 3-metric 
7ij, such that j'-^jjk = SI- 

There are many ways to formulate the Einstein equa- 
tions in an initial value, 3-1-1 form. The standard "ADM" 
3+1 formulation |l5j| writes the six space-space compo- 
nents of the Einstein equations, Eq. ^ as 12 equations 
that are first order in time 

^tltj = -2aKij + £^7y (5) 

LtK.j = a^''^R^j-2aK,''Kjk + aKK,j- (6) 

V{D,a + C^,K,j-a^^^R,j, (7) 

where Kij is the extrinsic curvature of a spacelike hy- 
persurface. Here, C is the Lie derivative operator, T^i is 
the covariant derivative operator compatible with the 3- 
metric 7^ , K is the trace of the extrinsic curvature, '■^•'-Rij 
is the 3-Ricci tensor, while ^'^'^Rij are the components of 
the 4-Ricci tensor that represent the perfect fluid source 
terms to the Einstein equations. The remaining four Ein- 
stein equations are the constraint equations, which are 
analytically satisfled on each of the spacelike hypersur- 
faces as long as they are satisfied on the initial slice. The 
Hamiltonian and momentum constraints are 

+ ^2 _ KijK'^ ~ 2a2G" (8) 

VjKh-'D.K - aG^^ = Q (9) 
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While one could base a numerical evolution code on the 
ADM equations, recent results (both empirical 0, IT^ 
and analytical studies 0,0, llo]) indicate that a more 
suitable choice would be the so-called BSSN formula- 
tion (16i. i17l i21i | . One feature of this formulation of the 
Einstein equations is that the 3-metric is decomposed 
into a conformal factor (j) and a conformal 3-metric 7^- as 



(10) 



where the determinant of the conformal 3-metric 7^ is 
identically 1 . Similarly, the extrinsic curvature is decom- 
posed into its trace and trace-free parts 



(11) 



where Aij is referred to as the conformal trace- free extrin- 
sic curvature, such that Aij^^^ = 0. In addition to the 
decomposition of the traditional ADM variables, a key 
ingredient in the BSSN formulation is the introduction 
of three new evolved variables, namely the three confor- 
mal connection functions F* 



(12) 



The final form of the evolution equations which we use 
in the numerical evolution of the Einstein field equations 
is given as 



1 



1 



dt 
dK 

'dt 

dt 
dt 



6 6 



(13) 

(14) 
(15) 



ae 

aKA 

^^,,VkV^a + e-^^C0{e^^A,,) 



3" 7y < 



'dt 



(16) 



1 



2aV]kA:^'' 



B. general relativistic hydrodynamics 



(17) 



We also perform a 3-)-l decomposition of the hydrody- 
namics equations, Eq. |31 Note that the 4- velocity is 



normalized u'^u^ = —1, so that its components can be 
written in terms of the three spatial velocity components 
as 



a 



(18) 



where W is the Lorentz factor W — \/ y/\ ~ jijv'-vL The 
specific enthalpy, /i, is given as 



h = l 



P/P, 



(19) 



where e is the specific internal energy density. 

The general relativistic hydrodynamics equations, 
Eq. 13 can be written in first order, flux conservative 
form as 



dtU + diF' = 5, 



(20) 



where the conservative hydrodynamical variables U are 
written in terms of the primitive variables {p, w',e} as 



^phW^Vj 
^{phW^ ~P- Wp) 

The flux vector F"^ is written as 





' D ' 










T 





(21) 



F' 



a (v' - /3'/a) D 



(22) 



and the source vector S is written as 







(23) 



C. discretization techniques 

We discretize each of the 3 spatial coordinate variables 
{a;, y, z\ using a constant spacing {Ax, Ay, Az}, e.g.. 



Xi = xq + i Ax, i = 0, • • • , ) 
We discretize the time coordinate t as 

tn+l = tn + At 



1. 



(24) 



(25) 



where we set At = 0.25 Ax for all dynamical simulations 
performed in this paper. 

Due to the fundamental differences in the phenomena 
being described by the Einstein field equations, Eq. ^ 
and the relativistic hydrodynamics equations, Eq. 13 the 
discretization methods that we employ for the two sets 
of equations are very different. In the case of the Ein- 
stein field equations, we expect the dynamical degrees of 
freedom to remain smooth and continuous. In the case 
of the relativistic hydrodynamical equations, we know 
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that shocks (discontinuities) can easily form in the phys- 
ical degrees of freedom. Thus, the discretization method 
used for the hydrodynamical equations will be more com- 
plicated in order to allow for the accurate treatment of 
shock propagation. The approach we use will be based 
on a finite differencing scheme employing High Resolu- 
tion Shock Capturing (HRSC) techniques. In order to 
use these techniques, a complete knowledge of the char- 
acteristic information is needed. We therefore require the 
eigenstructure of the Jacobian matrices in Eq.l^Ol namely 
dF'/dU for the flux vector defined in Eq.^ This is 
not a straightforward task, since the flux is expressed 
as a function of both the primitive and evolved hydrody- 
namical variables. What we require are a complete set of 
eigenvectors [fi] and corresponding eigenvalues Xi such 
that 



dF^ 



dU 



■,5. 



(26) 



(here, we present the spectral decomposition for the x- 
component of the Jacobian, since the decomposition for 
the other two spatial components of the Jacobian can be 
obtained by a straightforward permutation of the spa- 
tial coordinates {x,?/,z}). It turns out that the spectral 
decomposition contains a triply degenerate eigenvalue 



Ai = A2 = A3 = av^ - 



(27) 



A set of linearly independent vectors that span this de- 
generate space is given by 



ri 



hW{K-pcs^) 



,V,j.,Vy,V^, 1 



hW{K- pcs^) 



(28) 



f2 = [WVy, h{j^y + 2W'^V^Vy),h{jyy + 2W^VyVy), h{-fy, + 2W\yV,), VyW{2Wh " l)] , 

fg = [Wv,,h{^,, + 2W\^v,), h{jy, + 2W\v,),h{^,, + 2W\,v,),v,W{2Wh - 1)]^. 
The other two eigenvalues are given by 



A4 



1 — W^Cs^ 

with corresponding eigenvectors 



1, hW[v^- — ^— . , I , hWvy, hWv^, — ^-T. 1 



(29) 
(30) 

(31) 
(32) 



where the relativistic speed of sound in the fluid Cg is 
given by (see, e.g., 12^) 



dP 
dE 



X 
h 



P K 



P 



(33) 



We have set x 



dP 

dp 



and K 



dP I 

de I 



S is the entropy 



per particle and E is the total rest energy density which 
in our case is -B = p + pe. We use the above characteristic 
information to calculate the numerical fluxes {f*Y using 
the Piecewise-Parabolic Method (PPM), described in [23, 
01 . While the PPM method has been extended to special 
relativistic applications (see, e.g., H^), this is the first 
fully general relativistic application of the method. The 



discretization of the flux terms in Eq. 1201 are written 



(/*) 



dx 



i+l/2 



(/*) 



-1/2 



Aa 



0(Aa;2 



(34) 



In order to update the discretized hydrodynamical vari- 
ables, we simply perform a two-step predictor-corrector 
method, in order that the entire hydrodynamical update 
is done in a fully second order manner in both space and 
time (modulo the points where the hydrodynamical vari- 
ables obtain local extrema, where the accuracy of the spa- 
tial derivatives drop down to first order in space. This is 
a well known property of the so-called Godunov schemes; 
see, e.g., 

As previously stated, since the fields describing the 
gravitational degrees of freedom are expected to remain 
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smooth, we simply perform centered-in-space discretiz- 
ing of the spatial derivatives in Eqs. ^1 For dis- 
crete time evolution, we use the Iterated Crank- Nicholson 
method In order to achieve a completely second or- 
der method in both space and time for the coupled system 
of equations (the Einstein field equations and the rela- 
tivistic hydrodynamics equation), we use the time step- 
ping method described in Figure 

D. gauge choices and boundary conditions 

In the 3-1-1 initial value formulation of general relativ- 
ity, one is free to specify the slicing and spatial coordinate 
conditions by specifying the lapse a and shift respec- 
tively. The code described in the previous subsections 
has been written allowing for an arbitrary choice of these 
gauge variables. As described in the next section, each 
configuration in the CFQE-sequence approximation has 
a vanishing trace of the extrinsic curvature K . As we 
would like to compare our full general relativistic simula- 
tions to the CFQE-sequence approximation in an invari- 
ant manner, it is desirable to use the same slicing condi- 
tion during the full numerical simulation as that of the 
CFQE-sequence approximation. To this end, we would 
like to select the lapse function a such that the trace of 
the extrinsic curvature K remains 0. Notice that if one 
sets dK/dt = in Eq. then the equation becomes 
an elliptic equation for the lapse function a. We have 
thus implemented a multigrid solver (2&] for efficiently 
solving this equation, also known as the maximal slicing 
condition |29j |. However, solving an elliptic equation at 
every timestep can be numerically expensive. We there- 
fore also implement a variant of the so-called "l-flog" 
slicing condition for the lapse, 

^ = -2aK. (35) 

Note that this is a completely local condition, and is 
therefore computationally inexpensive. We use both the 
maximal slicing condition and the "l-|-log" slicing con- 
dition for simulations presented in this paper. For each 
result, we indicate which slicing condition is used. 

For the conditions on the shift, we use a slight mod- 
ification of the "Gamma- Freezing" shift equation [soj . 
Specifically, we implement the first integral form of the 
hyperbolic Gamma-driver (Eq. 46 of reference |3C|). 

^-Cir'-C2/?\ (36) 

where we set the constants Ci — C2 — 0.8 for all numer- 
ical simulations in this paper. We find that the Gamma- 
Freezing condition in this modified form, despite its sim- 
plicity, results in enhanced stability for the type of prob- 
lems studied in this paper, and was also very computa- 
tionally efhcient. 

The problem of boundary conditions in numerical rel- 
ativity is only now gaining the attention it deserves |3ll 




FIG. 1: A representation of tlie coupling between the hydro- 
dynamic predictor-corrector scheme (circles) and the iterated 
Crank-Nicholson method used for the integration of the Ein- 
stein field equations (squares). STEP 1: Simultaneous up- 
date of the general relativity and hydrodynamic equations via 
a Euler-predictor step (first order in time) to the half timestep 
n-l-1/2. STEP 2 through M: Update of the general relativ- 
ity equations via an iterative Crank-Nicholson scheme (second 
order accurate in time) to the n timestep, then compute 
a corrected n-\-l/2 state by averaging the n-\-l and n states. 
STEP M+1: Simultaneous update of the general relativity 
equations via a leapfrog step (second order in time) based on 
the n and n+1/2 states, and the hydrodynamics equations via 
the second half of the Euler-predictor step (first half applied 
in step 1) using a method of lines. STEP M+2: Update of 
the hydrodynamic equations to a virtual n 4- 2 timestep via 
a (first order in time) Euler-corrector step using method of 
lines. STEP M+3: A second order (in time) hydrodynam- 
ics update is obtained by averaging the corrected quantities 
of step M -I- 2 at level n-\-2 and the original variables at level 
n. 



[33, IS 0, nil ■ Here, we implement a simple, yet empiri- 
cally stable, boundary condition for the fields describing 
the gravitational degrees of freedom, Eqs. El due 
to Alcubierre |33|. Let f{t,x,y,z) represent a field on 
which we wish to impose this boundary condition in, say, 
the x-face of our computational boundary (the y- and z- 
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face are similar; averaging this method yields boundary 
conditions for the edge and corner points of the compu- 
tational boundary). We take as an ansatz the form for / 
as r — > cx) 

u{r — t) w{r + t) 



f{t,x,y,z) -> /oo + 



(37) 



Taking the derivative of this expression with respect to t 
and X, and eliminating u and u' yields 



^5/ _^ a/ 
r dt dx 



(38) 



where the function H — 2xw' . Eq.|2Hlis finite differenced 
in a second order fashion to obtain a boundary condition 
on the field /, where we have interpolated the function H 
from interior points, assuming a falloff of ^ for H . We 
adopt this boundary condition for all of the fields describ- 
ing the gravitational degrees of freedom. The bound- 
ary condition used for the hydrodynamical variables is a 
simple outflow boundary condition. Note that all of the 
hydrodynamical fields are trivially small at the compu- 
tational boundary (the atmosphere has a baryonic mass 
density that is 10^ times smaller than the central bary- 
onic mass density of the neutron stars), thus the bound- 
ary conditions used for these fields are relativity unim- 
portant. 



III. THE CONFORMALLY FLAT, 
QUASI-EQUILIBRIUM (CFQE) 
APPROXIMATION 

The mathematical assumptions that go into the con- 
formally flat, quasiequilibrium approximation to general 
relativity for binary, corotating neutron stars are as fol- 
lows (for the physical motivation behind the assumptions, 
see, e.g., 0). 

• The physical 3-metric 7ij is assumed to be confor- 
mally flat 

(39) 



7ij = il^'^S., 



I] ■ 



• The Lie derivative of the conformal metric '4'~'^lij 
with respect to the time variable t is identically 
zero: 



(40) 



• The trace of the extrinsic curvature and it's time 
derivative vanishes: 



K = 
Ct{K) - 0. 



(41) 
(42) 



• There exists an approximate timelike helical Killing 
vector field 



(43) 



• The 4- velocity of the fluid is proportional to the 
approximate Killing vector field <^'': 

(44) 

From Eq. I4UI the extrinsic curvature Kij takes the form 
1 /„ „ 2 



= ^ + ^J/^^ - ) (45) 

Using this, along with the assumptions of the conformally 
flat, quasiequilibrium approximation above, the Hamilto- 
nian constraint, Eq. |S1 and the momentum constraints, 
Eqs.El can be written as 



d'-d.i) + —K,jK'^ + 27riP^phW^ - P) = (46) 



and 



2a 



. 

respectively, where we define the conformal extrinsic cur- 
vature Kij as 



(48) 



Using this form of the Hamiltonian constraint, the max- 
imal slicing condition, CtK = 0, can be written as 

7a ~ ~ 

d'd^{atP)-—K^jK''+2TTi;*{2p{l+e)~3P-3phW^) = 0. 

(49) 

Since the quasiequilibrium approximation assumes the 
existence of a timelike Killing vector, we can analytically 
find the first integral of the relativistic Bernoulli equa- 
tion, which is simply 



— = const. 
h 



(50) 



for some constant fl. 



Along with the normalization condition for the 4- velocity, 
u'^Ufi = — 1, this equation can be explicitly written as 

(a^ _ ^4((^x _ ^^^2 ^ ^^y ^ ^^^2 ^ ^^2^2^ ^ ^^^^^ 

(51) 



A. CFQE configurations 

We use the algorithm detailed in 2] , along with a paral- 
lel multigrid solver, to simultaneously solve the algebraic 
Bernoulli integral equation, Eq. 1511 and the five elliptic 
equations corresponding to the Hamiltonian constraint 
fEa. l46|l . the three momentum constraints (Eas. l47|l . and 
the maximal slicing equation (Eg. I49|l . for the conformal 
factor tfj, the three components of the shift vector /?% the 
lapse a, and the matter fields, thus producing a single 
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FIG. 2: The ADM mass and baryonic mass (in units of Mq) 
as a function of central density (in units of nuclear density) for 
single, static neutron star configurations. All studies in this 
paper are done with stars of baryonic mass 1A9Mq (square). 
The corresponding ADM mass of this single static configura- 
tion is IAMq (circle). 



CFQE configuration. In numerically solving these ellip- 
tic equations, we use identical boundary conditions to 
those in 0^ . We have also assumed a polytropic equation 
of state, 



P = (r - l)pe = kp^ 



(52) 



where k is the polytropic constant and T is the adiabatic 
index. The baryonic mass of each neutron star in a CFQE 
configuration, Mq, is defined as half of the total rest mass 
of the configuration. 



1 



(fx ^pW. 



(53) 



In this paper, we set F = 2 and k = 0.0445|-, where p„ 
is nuclear density (approximately 2.3 x 10^^ g/cm?). For 
these values of parameters, a single static neutron star 
configuration that is stable has a maximum ADM mass 
of 1.79Mq and a baryonic mass of l.OTM© {Mq is 1 solar 
mass). For the studies in this paper, we use neutron stars 
with a baryonic mass of 1.49Mq each, which is approx- 
imately 75% that of the maximum stable configuration. 
The ADM mass of a single static neutron star for this 
configuration is 1.4M0 (see Figure|2l. 

Once the baryonic mass of each of the stars is fixed, 
the only remaining degree of freedom in the specification 
of a CFQE corotating configuration is the separation of 
the two stars. A natural invariant way of specifying the 
separation of the two neutron stars is to calculate the 
geodesic distance on the constant t hypersurface between 
the points in each of the neutron stars that corresponds 
to the maximum baryonic mass density. Specifically, if 



and 



'M2 



are the spatial coordinates of the points 



of maximum baryonic mass density in the first and sec- 
ond star, respectively, then the geodesic whose length 



will represent the (spatially) invariant separation of the 
neutron stars is the curve 



where 



X'(A), Ae[0,l], 



a:*(a=o) = 4^1, x\\=\) = x 



M2; 



such that 



F' 



dX^ dX^ 
dX dX 



0. 



(54) 



(55) 



(56) 



The geodesic distance £1^2 between the two stars is now 
defined as 



\=i 



dX 



A=0 



dX^ dX3 



dX dX 



(57) 



We use a relaxation technique for numerically solving 
the geodesic ODE, Eq. [SHI and then use a standard 
quadrature formula for evaluating the integral for ^1^2, 
Eq. 123 We choose the discretization for solving Eq. l5?il 
and Eq.l^to be 10 times as fine as the 3D discretization, 
and use a third order accurate interpolator to calculate 
the Christoffel symbols F'jfe from the 3D grid. The max- 
imum baryonic mass points, x\j-^ and x\,j2^ are located 
by finding the maximum of a 3D third order accurate in- 
terpolation polynomial, the interpolation being centered 
on the discrete point on the 3D grid that has the largest 
value of the baryonic mass density in each star. 

In Figure|31 we plot logarithmically spaced contours of 
the baryonic mass density in the equatorial (z = 0) plane 
for two CFQE configurations representing the smallest 
and largest separations used as initial data for the 3D 
numerical simulations performed in this paper. 



B. The CFQE-sequence approximation 

A CFQE-sequence is constructed by stringing together 
several constant baryonic mass CFQE configurations, 
each configuration differing in only the separation of the 
neutron stars. The idea is that, due to the fact that 
the time scale of the gravitational radiation, and thus, 
that of the orbital decay, is much longer than the orbital 
timescales, the binary neutron stars are considered to 
be in "quasiequilibrium" . Thus, at any particular time, 
a binary neutron star system can be described by one 
particular CFQE configuration; the effect of the gravi- 
tational radiation is, over the timescales of one orbit, to 
alter the configuration to a new CFQE configuration with 
a slightly smaller separation. 

It is typical in CFQE approximation studies to cal- 
culate an effective binding energy for each configura- 
tion. Following Q, we define the (dimensionless) effective 
binding energy Ef, of a single configuration to be 
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Madm - 2M_ 
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(58) 
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FIG. 3: Contours of the baryonic mass density for two CFQE 
configurations, taken in the equatorial {z — 0) plane. The 
ADM mass of each neutron star in isolation is IAMq. The 
value of the mass density for the first contour is 0.9 that of the 
maximum mass density of the star, with the others decreasing 
by a factor of 0.5 each. The CFQE configuration in the top 
panel has a geodesic separation of = 23.44 (which is 

51.70 km), and the configuration in the bottom panel has a 
geodesic separation of = 35.72 (which is 79.00 km). 

These configurations represent the smallest and largest sepa- 
ration CFQE configurations that are used as initial data for 
the fully consistent general relativistic numerical calculations 
done in this paper. 

where Madm is the ADM mass of the configuration, and 
Mmsoq is the ADM mass of a single static neutron star 
in isolation with rest mass Mq. The usual expression for 
the ADM mass of an asymptotically fiat spatial slice is 

Madm = ^ lim I dA {^-^ - N\ 

(59) 

where iV* is the unit outward normal to the sphere of con- 
stant radius r, which is the domain of integration. Using 
the fact that the CFQE configurations are conformally 
flat, the ADM mass of any configuration reduces to the 
volume integral 

Madm = j d\ {d'd.'if) . (60) 

In Figure 0] we plot the effective binding energy Ei, 
as a function of the angular velocity parameter for a 
constant baryonic mass CFQE-sequence (each star has a 
baryonic mass of I.49M0). The highest value of for 
each of our calculated CFQE-sequences corresponds to a 
CFQE configuration where the neutron stars are touch- 
ing. The configuration corresponding to the minimum 
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FIG. 4: The effective binding energy iSj, as a function of the 
orbital angular velocity parameter Q. for constant baryonic 
mass CFQE-sequences. The baryonic mass of each of the 
neutron stars is 1.49M0. Shown are results for our CFQE 
configuration solver at different resolutions and outer bound- 
ary placements. Tables IV, V, and VI of reference 2j were 
interpolated to obtain the results from Baumgarte, et. al. 

effective binding energy E), is the ISCO configuration. 
It is this point along the CFQE-sequence approxima- 
tion where the quasiequilibrium configuration goes sec- 
ularly unstable; the subsequent evolution of the system 
is thought to enter a "plunge phase" where the neutron 
stars coalesce within an orbital timescale. 

Of course, when using a computer code to solve differ- 
ential equations, results produced at one single resolution 
are meaningless due to the fact that we have no way of 
knowing how big the numerical errors are. It is always 
important to make numerical calculations at different res- 
olutions, in order to be able to assess the numerical errors 
that are inherent to any discretization. Not only do we 
have the usual truncation errors (which are due to the 
fact that higher order terms in the Taylor expansion of 
functions have been dropped in our finite difference ap- 
proximation) of any finite difference approximation, we 
also have boundary errors, as we are solving a set of 5 
coupled elliptic equations to arrive at any CFQE configu- 
ration. Therefore, we must make at least three numerical 
experiments in order to assess the effect of both resolu- 
tion and boundary placement on the numerical results. 
We show our CFQE configuration results for three com- 
binations of resolution and outer boundary placement, 
alongside the results obtained by Baumgarte, et. al. 0]. 
Our CFQE sequence results are in fairly good agreement 
with those of P| , considering the two terms that are sub- 
tracted in the construction of Eb, Eq. |2H1 are the same 
to the first two or three significant digits. The resolu- 
tion used in the two studies arc comparable. However, 
in , the number of gridpoints across the neutron star 
was kept fixed. We, on the other hand, keep (for each 
sequence calculated) the boundary fixed, and thus have 
different numbers of gridpoints across the star in each 
sequence. Notice that, as we both increase the resolution 



and increase the distance from the center of mass to the 
boundary of the computational domain, the shape of the 
curves in Figure 01 remain roughly the same. In fact, we 
can perform a Richardson extrapolation on our numeri- 
cal results by fitting each point on the curve to an error 
function. For example, we can use the error function 



'numerical 



i = (-B6)exact + Ci(Ax)' + -%, (61) 



where {Eb)numericai IS the numerical value for the bind- 
ing energy Ei, obtained using a specific discretization Ax 
and coordinate distance from the center of mass to the 
outer boundary rid (we use rid to denote this distance 
in solving the elliptic equations for the CFQE configura- 
tions, and we will use rb to denote the coordinate distance 
from the center of mass to the outer boundary of the 
computational domain used during the dynamical evolu- 
tions. Of course, we always have rid > fb-). Note that 
the first nonzero term for the expansion of the bound- 
ary error does not contain a l/rid term. This is due to 
that fact that the boundary conditions we use in solving 
the elliptic equations for the CFQE configurations, which 
are identical to those used in are exact to this order. 
{Eb) exact is the binding energy in the limit as Ax and 
rid — > oo, i.e. that given by an exact solution to the dif- 
ferential equations. Of course, it is {Eb)exact that we are 
interested in. We use our three CFQE-sequences shown 
in Figure ^ to solve for the three unknowns {Eb)exact, 
Ci , and C2 from Eq. |^ A generous estimate of the to- 
tal error, which will be used for the size of the error bars, 
is 



error = max{ | Ci (Ax)^^^^^^ | , 



The results are plotted in Figure |S| 



(62) 



C. neutron star spin corrections to the 
CFQE-sequence approximation 

Note in Figure El that the effective binding energy at- 
tains a minimum as the separation between the neutron 
stars decreases. As explained in the introduction, the 
CFQE-sequence approximation thus predicts a secular 
instability at this orbital separation, at which point the 
evolution of the system will change from quasistationary 
into a plunge phase, where the neutron stars would merge 
on timescales of the orbital period. Several types of argu- 
ments have been used to support this claim ^^]. On the 
one hand, it is intuitively clear that, in the presence of 
a dissipation mechanism (gravitational waves are slowly 
dissipating the binding energy of the binary system), if 
the binding energy would increase with decreasing sep- 
aration, this would energetically be an unstable situa- 
tion. There are also turning point theorems 0, 113. IT^ in 
full nonlinear general relativity which state that if a se- 
quence of equilibrium configurations attain a minimum 
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FIG. 5: The effective binding energy Ef^ as a function of the 
orbital angular velocity parameter fl. The data in this plot 
were obtained by using Richardson extrapolation from the 
data in Figure |11 The form of the truncation and boundary 
error are given by Eg. 1611 The error bars were computed using 
Eg. 1621 These results represent what one would obtain in the 
limit as both Aa; and the location of the outer boundary 
goes to 00. 



in the ADM mass-energy, then this equilibrium config- 
uration is secularly unstable. Hence, to the accuracy of 
the CFQE-sequence approximation, the minimum in Fig- 
ure [3 is often referred to as the innermost stable circular 
orbit (ISCO) configuration. 

One problem with the view that these CFQE- 
sequences somehow approximate fully consistent solu- 
tions to the Einstein equations coupled to the relativistic 
hydrodynamic equations has to do with the spin of the 
the individual neutron stars. Note that in each of the 
CFQE configurations, the 4- velocity of the fiuid is as- 
sumed proportional to the timelike helical Killing vector 
field. The neutron stars are thus spinning at the exact 
same frequency as the orbital frequency, Vl. The neutron 
stars are said to be corotating with the orbital motion. 
However, it has been known for some time that realistic 
binary neutron stars cannot be tidally locked during the 
late stages of inspiral (i.e., less than 1000 orbits until the 
final plunge) [36j- 

This raises a question: if a CFQE configuration is used 
as initial data, should we expect the subsequent solution 
to the Einstein equations to follow the CFQE-sequence 
approximation, and keep the stars tidally locked as the 
stars slowly inspiral? In light of the results of [s^, the 
answer must be no. To first order, we would expect the 
stars to retain approximately the same spin during the 
dynamical evolutions as is given in the CFQE configura- 
tion that was used as initial data. 

The resolution of this first question presents a second: 
what is the status of the ISCO? If the stars do not, in fact, 
stay tidally locked during dynamical evolution, will this 
affect the location and/or existence of a turning point in 
the binding energy curve, i.e. the CFQE ISCO configu- 
ration? One can gain insight into this question by com- 
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FIG. 6: The effective binding energy E'l,, Eq.|S31 as a function 
of the orbital angular velocity parameter Q. The method for 
obtaining the data and error bars in this plot are equivalent 
to those used in producing Figure |K1 



paring the spin kinetic energy of relativistic, uniformly 
rotating neutron stars to the energy scales involved in 
calculating the binding energy, Eq. 1581 To this end, we 
define a new effective binding energy that explicitly takes 
the spin kinetic energy of the neutron stars into account: 



MadM - ^MnSoo - 2 AMrnSoo{^) 



where AMjinsoo{^) is the difference between the ADM 
mass of an isolated stationary neutron star with baryonic 
mass Mo uniformly rotating with angular velocity il and 
the ADM mass of an isolated static non-rotating neu- 
tron star with baryonic mass Mq- Thus, AMfi^soci^) 
represents the spin kinetic energy of a uniformly rotating 
neutron star rotating with an angular velocity f2. We plot 
this new effective binding energy i?^ in Figure El As can 
be seen, the binding energy where the spin kinetic energy 
of each of the neutron stars has been manually factored 
out no longer attains a minimum (the data point with the 
highest value of corresponds to the CFQE configura- 
tion in which the two stars are touching) . This indicates 
that the secular stability of very close binary neutron 
stars could depend very sensitively on the spin-orbit cou- 
pling of the binary system. The two extreme cases, that 
of tidally locked neutron stars and non-rotating neutron 
stars, are depicted by Figures |31 and El respectively. The 
fact that we expect a weak spin-orbit coupling suggests 
that, in the fully consistent general relativistic hydrody- 
namics simulations, where we do not expect the neutron 
stars to remain tidally locked, there may be no turning 
point and that there may not, in fact, be any ISCO con- 
figuration. 



IV. COMPARING NUMERICAL EVOLUTION 
IN FULL GENERAL RELATIVISTIC THEORY 
WITH THE CFQE-SEQUENCE 
APPROXIMATION 

Using our numerical evolution code described in Sec- 
tional along with the CFQE configurations described in 
Section UTTI as initial data, we analyze to what degree the 
CFQE configurations actually correspond to quasiequi- 
librium configurations in full general relativity. We do 
this by computing fully consistent general relativistic evo- 
lutions numerically using CFQE configurations as initial 
data, and comparing the resulting solution of the Einstein 
field equations with the spacetime obtained by the CFQE 
approximation. Recall that one major assumption in the 
CFQE approximation method is that the change of the 
binary separation is small on timescales less than one or- 
bital period. If the solution to the Einstein equations ob- 
tained by using any particular CFQE configuration as ini- 
tial data deviates significantly from the CFQE-sequence 
spacetime on suborbital timescales, then that particu- 
lar CFQE configuration cannot be thought of as repre- 
senting a quasiequilibrium configuration. Of course, the 
quasiequilibrium approximation becomes better as the 
separation of the neutron stars increases. That question 
then becomes: how close can the neutron stars be before 
the CFQE approximation breaks down? In this section, 
we describe a generic method that we subsequently use to 
answer this question. In short, we use CFQE configura- 
tions of differing initial separation as initial data for our 
fully consistent general relativistic code, and compare the 
resulting solutions to the Einstein field equations to the 
CFQE-sequence spacetime in a coordinate independent 
fashion. The magnitude of these differences, which we 
find to decrease as the initial separation of the neutron 
stars increases (as expected), will therefore quantify the 
degree of failure of the CFQE configurations to be truly 
quasiequilibrium configurations in full general relativity. 

Recall that the CFQE-sequence approximation is done 
in a maximally sliced {K = 0) coordinate system. In 
order to facilitate a meaningful comparison between 
the fully general relativistic calculation and the CFQE- 
sequence approach, we adopt this slicing condition for all 
numerical computations done in this section. Thus, we 
determine the lapse function a by solving the maximal 
slicing condition, which is the elliptic equation obtained 
by setting dK/dt to in Eq.^^ In order to decrease the 
computational cost involved in solving this elliptic equa- 
tion, we only solve the maximal slicing condition every 
eight timesteps. We use "1 -I- log" evolution. Eg. 1351 for 
the lapse a for the timesteps where we do not solve the 
maximal slicing condition. We have performed several 
numerical tests where we do, in fact, solve the maximal 
slicing condition for every timestep, and have seen only 
a negligible affect on the results (in fact, the measures of 
the violation of the Einstein equations by the CFQE ap- 
proximation that we find in this section slightly increase 
when we solve the maximal slicing condition for the lapse 
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Config. 


Mo/Mq 


^1,2/A/o 


J /Mo' 


Madm/Mo 




NS-1 


1.490 


23.44 


3.770 


1.857 


0.01547 


NS-2 


1.490 


25.94 


3.762 


1.857 


0.01296 


NS-3 


1.490 


29.78 


3.812 


1.858 


0.01022 


NS-4 


1.490 


35.72 


3.934 


1.859 


0.007460 



TABLE L CFQE configuration parameters used as initial 
data for fully consistent general relativistic simulations in Sec- 
tion ^3 For reference, the ISCO configuration has a proper 
geodesic separation of ^1,2 = 24.0 Mq. 



a at every timestep). 

In order to meaningfully compare our fully consistent 
general relativistic simulations with the CFQE-sequence 
approximation, we must compare quantities that are in- 
dependent of our choice of coordinates. Since we are us- 
ing the same slicing condition as that assumed by the 
CFQE-sequence approximation (namely, maximal slic- 
ing), we must compare our numerical results with the 
CFQE-sequence approximation using quantities that are 
either fully 4-invariant quantities, or 3-invariant quan- 
tities (quantities that are invariant under spatial co- 
ordinate transformations). In the following, we will 
construct measures for comparing our numerically con- 
structed spacetime with the CFQE-sequence spacetime 
based on both the conformal flatness assumption and the 
assumption of the existence of a timelike helical Killing 
vector field. Both of these measures will be exactly zero 
for the CFQE-sequence spacetime, but will not neces- 
sarily be zero for the full solution to the Einstein field 
equations using CFQE configurations as initial data. The 
magnitude of these measures will thus quantify the mag- 
nitude of the CFQE configuration's failure to correspond 
to a true quasiequilibrium configuration in full general 
relativity. In Table ^ we show the parameters for the 
CFQE configurations that are used as initial data for the 
simulations done in this section. 



A. The Killing vector field assumption of the 
CFQE-sequence approximation 

Recall from Section IIIII that one assumption of the 
CFQE-sequence approximation is the existence of an ap- 
proximate timelike helical Killing vector field, Eq.^Jl As 
we are assuming a corotating matter field, the 4- velocity 
of the fluid must be proportional to this Killing vector 
field, Eq. ^2 Here, we monitor this quasi-equilibrium 
(QE) assumption in our fully consistent general relativis- 
tic numerical calculations. 

That the 4- velocity u°' be proportional to a Killing vec- 
tor field is equivalent to the vanishing of a symmetric, 
type-(0,2) 4-tensor Qat 

Qab = ^aUb + VbUa + UaUb + UbUa (64) 

where a° = u^Wbu"" is the 4-acceleration of the fluid 



(Va denotes the covariant derivative operator compat- 
ible with the 4-metric gab)- Notice that the quantity 
Qabu"" vanishes identically. We can thus monitor the 
space-space components of Qab during our simulations 
as a way of monitoring how well the 4-velocity stays 
proportional to a Killing vector field. Define Qij as the 
projection of Qab onto the constant t spatial slice: 



?ij — Pi Pj Qab — Qab 



(65) 



where Fab = gab + riarib is the projection operator onto 
the constant t spatial slices. The unit normal to these 
spatial slices is = at^ + In our Cartesian coor- 
dinates a;* for the spatial slices, the components of Qij 
form a 3x3 matrix. The norm of this matrix, which it- 
self is a coordinate-independent quantity, is the square 
root of the largest eigenvalue of QijQ\, where we have 
raised and lowered 3-indices with the 3-metric. In our 
case, Qij is symmetric, and the matrix norm reduces 
to the largest eigenvalue of Qij itself. We denote this 
coordinate-independent value of the norm of Qij as \Qij\. 
Note that if the 4-velocity is proportional to an exact 
Killing vector, \Qij\ will be exactly zero. Of course, in 
the fully consistent general relativistic treatment, \Qij\ 
will not vanish. What is required is a sense of the relative 
size of \Qij\- Notice that in Eq. |^ Qab is constructed 
out of only two separate (symmetric) parts, the Vm part 
and the ua part. We can therefore naturally normalize 
\Qij\ by the norms of these two principle parts. If we 
define 



Qlab = ^aUb + VfcUa 
Q2ab = Uaab + Ubtta, 



(66) 
(67) 



then a naturally normalized scalar field Q which denotes 
the deviation from the 4-velocity being proportional 
to a Killing vector field is 



\Q^. 



max{ Qij- , Q2y } 



(68) 



This normalization provides a measure for the deviation 
of the 4-velocity from being proportional to a Killing 
vector field (the QE assumption of the CFQE approach); 
a value of Q = signifies that the 4-velocity of the fiuid 
is exactly proportional to a timelike Killing vector, while 
a value of Q of order unity would signify a significant 
violation of the QE assumption. The monitoring of Q 
during a fully consistent general relativistic simulation is 
then a quantitative measure of the accuracy of the QE 
approximation. Since Q is meaningful only inside the 
fiuid bodies, a natural global measure of the magnitude 
of Q is its the baryonic mass weighted integral, denoted 
by (Q): 



(Q) 



Jd^x \Q\^pW 
J (Px ^pW ' 



(69) 
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FIG. 7: We plot the quantity (Q) (Eq. 169^ as a function of 
time for fully consistent general relativistic numerical calcu- 
lations, using CFQE configuration NS-1 as initial data (see 
Table HJ. A variety of discretization parameters A.x are used. 
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FIG. 8: We plot the quantity (Q) (Eq. as a function of 
time for fully consistent general relativistic numerical calcu- 
lations, using CFQE configuration NS-2 as initial data (see 
Table HJ. A variety of discretization parameters Aa; are used. 



where the integrals are taken to be over the entire spatial 
slice. 

In Figure [TJ we plot (Q) as a function of time for a 
small fraction of an orbit in the fully consistent general 
relativistic numerical simulations. Configuration NS-1 
was used as initial data (see Table Pi . Various resolu- 
tions were used, along with different numbers of grid- 
points for the computational domain. As stated in the 
section where we numerically solved for the CFQE con- 
figurations. Section [ill Bl it is important to run any sim- 
ulation at multiple resolutions and boundary placements, 
in order to assess the magnitude of the boundary error 
and the finite difference truncation error on the numerical 
results. 

Notice in Figure [7| that the value of (Q) appears to be 
converging to a curve that attains a maximum value of 
approximately {Q) = 0.26 after 1.5% of an orbit {2tt/VI 
is approximately 1 orbital period). This is quite a fast 
growth, and a (perhaps unexpectedly) large value at- 
tained in such a short time. Recall that (Q) is normalized 
such that a value of {Q) of order unity represents a sig- 
nificant failure of the 4- velocity to be proportional to 
a Killing vector. One reason for this rapid growth has 
been discussed in Section llll CI there is no mechanism in 
the actual evolution to spin up the neutron stars in order 
that they remain in corotation with the orbital motion. 
This corotation condition is "force-fed" into the CFQE 
configurations. Of course, one expects that the CFQE- 
sequence approximation should become better as the sep- 
aration between the two neutron stars increases. The or- 
bital angular velocity, and thus, the spin of the neutron 
stars, decreases with increasing separation. Hence, both 
the gravitational radiation reaction and the error intro- 
duced by the artificial spin-up of the individual stars due 
to the corotation assumption are lessened. These two 
factors will increase the validity of CFQE approximation 
for increased neutron star separations. 

In Figures |H1 El and ^| we plot {Q) as a function of 
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FIG. 9: We plot the quantity (Q) (Eq. IHIIjl as a function of 
time for fully consistent general relativistic numerical calcu- 
lations, using CFQE configuration NS-3 as initial data (see 
Table HJ. A variety of discretization parameters Ax are used. 



time for increasing neutron star binary separations. No- 
tice that, for each separation, the shape of the curves 
appear to be converging to something that is qualita- 
tively similar to that of Figure [7| a quick increase to a 
maximum value. Also note that this maximum value is, 
as expected, decreasing for increasing neutron star sepa- 
ration. The best resolution for the maximum initial neu- 
tron star separation (NS-4) in Figure ^| (the solid line) 
has (Q) attaining a maximum value that is already lower 
than 0.2. A natural question then arises: can one predict 
how far the initial separation of the neutron stars should 
be in order that the maximum value of (Q) obtained in a 
short time scale in the consistent general relativistic the- 
ory be bounded by some number, say, (Q) = 0.1? One 
problem is in the details of the numerical simulations: as 
the neutron star separation increases, the computational 
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FIG. 10: We plot the quantity (Q) (Eq. 1691 as a function 
of time for fully consistent general relativistic numerical cal- 
culations, using CFQE configuration NS-4 as initial data (see 
Table HJ. A variety of discretization parameters Ax are used. 




FIG. 11: The maximum value of (Q) (Ea. l69ll obtained in our 
fully consistent general relativistic simulations using different 
initial CFQE configurations with initial separations £i_2. The 
error bars are computed using the error function Eg. lTUI 



resources demanded by the problem become larger. In 
other words, for a given fixed amount of computational 
resources, the ability to resolve each neutron star (e.g., 
the number of discrete gridpoints across each star) is di- 
minished as the initial separation of the neutron stars is 
increased. This can be seen directly when comparing the 
resolutions used for the simulations performed for Fig- 
uresOIHlEl and llUI as the initial neutron star separation 
is increased in going from CFQE configuration NS-1 to 
NS-4, the resolution necessarily decreases. These effects 
can be seen directly in the size of the error bars of Fig- 
ure ^2 where the maximum value of (Q) obtained in the 
short term general relativistic simulations is plotted as a 
function of the initial geodesic separation ti^2- For each 
data point, we use the maximum value attained by the 
highest resolution curve (solid line) in each of Figures [7|- 
IIUI As always, it is important to estimate the truncation 
and boundary errors in any numerical calculation. Here, 
we use an error estimate of the form 

max) numerical max ) exact + 

Ci(Ax)+C2(Aa;)V-%.(70) 

n 

Note that we now include a term that is linear in the 
discretization parameter Ax. This is due to the fact 
that, although we are using second order methods for 
the discretization of the Einstein equations, the HRSC 
methods used for solving the relativistic hydrodynamics 
equations (described in Section III C|) are only first order 
(in space) accurate at points where the hydrodynamical 
variables obtain a local extrema. Also, while we expect 
the outer boundary condition on general dynamical sim- 
ulations to be better represented by a 1/rt, error term, we 
note that in these short time scale simulations, the neu- 
tron stars are not even causally connected to the outer 
boundary. The only boundary error in the calculation is 
that due to the initial data CFQE configuration solve. 



whose boundary error decreases as The error bars 

used in Figure 1111 are computed using the 4 numerical 
results obtained by varying the resolution and boundary 
placement for each configuration, and solving Eq. I7UI for 
the constants {{Q) max) exact, Ci, C2, and C3. The size 
of the error bar is then set equal to the largest of the 
absolute value of each individual error term in Eq. 1701 
With this generous estimate of the error of the maxi- 
mum value of {Q), we can put a lower bound on the ini- 
tial geodesic separation of the CFQE configuration that 
must be used as initial data in fully general relativistic 
simulations such that the Killing vector field assumption 
is valid, e.g., (Q) < 0.1. While the value of 0.1 may be 
somewhat arbitrary (one may desire an even more strin- 
gent criterion), the method we use to analyze the CFQE 
data is quite general. In Figure 1111 we fit a reciprocal 
power law, ^ + + to both the 

maximum value obtained in our fully consistent general 
relativistic calculations, as well as to the lower bound of 
the estimated error in our calculation. We can see that 
one would have to use a CFQE configuration initial data 
set with a geodesic separation between the neutron stars 
of at least £1.2 = 46.8 Mq in order for the subsequent 
solution to the full Einstein field equations to satisfy the 
Killing field assumption of the CFQE-sequence approxi- 
mation to 1 part in 10 (this separation parameter could 
actually be as low as ^i_2 = 41.2Mo, taking into account 
the errors of our calculations, see Figure ITT)) . This sep- 
aration corresponds to roughly twice that of the ISCO 
separation. 
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FIG. 12: We plot the quantity {H)p (Eq.^IJ as a function 
of time for fully consistent general relativistic numerical cal- 
culations, using CFQE configuration NS-1 as initial data (see 
Table HJ. A variety of discretization parameters Ax are used. 



B. The conformal flatness assumption of the 
CFQE-sequence approximation 

One other assumption in the CFQE-sequence approx- 
imation is that of conformal flatness (CF). It is often 
argued that this assumption is, in some sense, equiva- 
lent to assuming that there is no gravitational radiation 
in the configuration. However, the statement that con- 
formally flat configurations have zero gravitational radi- 
ation content is very questionable, especially in the case 
of the CFQE-sequence approximation which is not even 
consistent with the full set of Einstein equations. Here, 
we analyze this CF assumption in full general relativ- 
ity. We start with a CFQE configuration as initial data, 
and perform fully consistent general relativistic numeri- 
cal evolutions, monitoring the conformal flatness of the 
spatial slices. As we are using the same slicing condition 
(maximal slicing) as that of the CFQE-sequence approx- 
imation, we only require a 3-invariant that will allow us 
to monitor the conformal flatness assumption during the 
simulations in a coordinate independent way. The 3-Bach 
tensor is one such 3-invariant. It is defined on the spatial 
slice as 
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and can be shown to vanish if and only if the 3-metric 



7ij is conformally flat. The Cotton- York tensor, 
related to the 3-Bach tensor by 
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where eijk is the natural volume element 3-form. We de- 
fine the scalar H as the matrix norm of the Cotton- York 
tensor, normalized by the size of the covariant derivative 



of the 3-Ricci tensor: 
H 



\H„ 



(73) 



where, just as in the previous section, \Hij\ denotes the 
matrix norm of the components of Hij in our Cartesian 
coordinates. Note that H vanishes on conformally flat 
spatial slices, and is normalized to provide a local mea- 
sure for determining how much the spatial slice is de- 
viating from conformal flatness. For a global measure, 
we define the baryonic density weighted norm denoted as 



^ Jd^x \H\^pW 



(74) 



where the integrals are taken to be over the entire spatial 
slice. 

In Figure El we analyze the CF assumption of the 
CFQE-sequence approximation by plotting {H)p as a 
function of time for fully consistent general relativistic 
simulations. We use for initial data the CFQE configura- 
tion NS-1 (see TablePl. As it is always imperative to run 
a numerical code at multiple resolutions and boundary 
placements to assess the numerical errors, we use a vari- 
ety of discretizations and grid sizes. Using this measure 
of the violation of the conformal fiatness assumption in 
the CFQE-sequence approximation, we see that the as- 
sumption holds to roughly 1 part in 20, for this initial 
data. We can also see that this measure, as compared to 
the measure for the QE assumption (Eq.lHHI Figurc[7}, is 
not as sensitive to resolution. In other words, the numer- 
ical truncation error for this particular measure is not as 
large. 

Again, we would expect that the CF assumption to be 
better for larger initial neutron star separation. In Fig- 
ures ^1 El and 1 151 we plot the measure of the violation 
of the conformal flatness assumption {H)p, Eq. El in 
our general relativistic simulations using CFQE configu- 
rations NS-2, NS-3, and NS-4, respectively, as initial data 
(see Table nj. We can see that the violation of the con- 
formal fiatness assumption does, in fact, decrease with 
increasing initial neutron star separation. As with the 
QE assumption, we can use the results of Figures [T^- [TBI 
to predict the initial neutron star separation one would 
need for a CFQE configuration to satisfy the CF assump- 
tion to some prescribed tolerance. For example, we may 
want to start our general relativistic calculations with ini- 
tial data that corresponds to a CFQE configuration such 
that the error in the CF assumption is below one part in 
100, as measured by the quantity {H)p (Eg. I74|l . In Fig- 
ure El we plot the maximum value of quantity {H)p as 
a function of initial geodesic separation £1,2 attained in 
our fully consistent general relativistic numerical simula- 
tions using the four CFQE configurations from TableQlas 
initial data (see Figures El - El • Again, we use Eg. [701 
and the method described in Section IIV Al to compute 
the numerical errors (both truncation errors and bound- 
ary errors) made in the calculation. We fit an inverse 
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FIG. 13: We plot the quantity {H)p (Eq.^IJ as a function 
of time for fully consistent general relativistic numerical cal- 
culations, using CFQE configuration NS-2 as initial data (see 
Table HJ. A variety of discretization parameters Ax are used. 
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FIG. 15: We plot the quantity (Ji")p (Eq. |^ as a function 
of time for fully consistent general relativistic numerical cal- 
culations, using CFQE configuration NS-4 as initial data (see 
Table HJ. A variety of discretization parameters Aa; are used. 
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FIG. 14: We plot the quantity (iy)p (Eq. as a function 
of time for fully consistent general relativistic numerical cal- 
culations, using CFQE configuration NS-3 as initial data (see 
Table HJ. A variety of discretization parameters Ax are used. 



FIG. 16: The maximum value of {li)p (Eq. I74|l obtained 
in our fully consistent general relativistic simulations using 
different initial CFQE configurations with initial separations 
The error bars are computed using the error function 
Eg. 1701 



V. LONG-TERM GENERAL RELATIVISTIC 
NUMERICAL SIMULATIONS 



power series function -f -1- + rr^ to 

the four data points, as well as to the lower bound of the 
error, in Figure 1161 As can be seen, one would have to 
use initial data corresponding to a CFQE configuration 
with neutron star geodesic separation of approximately 
46.7Mo or greater in order for (i?)p to be 0.01 or less 
in the subsequent solution to the Einstein field equations 
coupled to the hydrodynamics equations. Recall from 
the previous section that this separation would also sat- 
isfy the Killing field assumption of the CFQE-sequence 
approximation to 10%. 



In the previous section, we performed many short 
timescale (sub-orbital) general relativistic simulations us- 
ing CFQE configurations as initial data. There, the focus 
was on determining the intrinsic error in using the CFQE 
configurations as initial data to model astrophysical neu- 
tron star binaries. Here, we study the long timescale 
(i.e., multiple orbital timescale) evolution of the system 
in order to investigate the stability and accuracy of long 
term fully general relativistic numerical integrations of 
the neutron star binary system. 

There are several reasons why it is extremely difficult 
to study large separation binary neutron stars on time 
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scales longer than the orbital period in full numerical 
relativity in an accurate fashion. One basic reason is ap- 
parent from the dynamics of two Newtonian point masses 
in circular motion. The orbital period Torb increases as 
the separation D^ep of the masses as 
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(75) 



Thus, any simulation of neutron stars with larger sepa- 
ration will naturally have to be run to longer times in 
order to capture a single orbital period. Also, from the 
standard quadrapole formula, the energy in gravitational 
waves Eg emitted per unit time for particles in circular 
motion decreases with increasing separation Dgep as 



dt ^^'P' ■ 



(76) 



Obviously, the time for any simulation to track the orbital 
motion of compact binaries through the plunge phase to 
the final merger will be extremely sensitive to the initial 
separation. 

We can get some idea of the computational resources 
needed to accurately simulate binary neutron stars with 
initial geodesic proper separations of between 25 Mq 
and 35 Mq (which corresponds to roughly 2.5 Rns and 
4.5 i?Ars, where Rns is the neutron star radius) from the 
results of the previous section. In Section IIVI we per- 
formed general relativistic numerical calculations using 
roughly these separations (see Table for only several 
percent of one orbital period. It took numerical configu- 
rations of over 500'^ gridpoints in order to obtain resolu- 
tions high enough for a confident prediction of the error 
of the simulations. In order to perform simulations on or- 
bital timescales, we need to increase the simulation times 
by two orders of magnitude. While this is already quite 
difficult, the situation is even more demanding: such a 
simulation time is much greater than the light crossing 
time of our computational domain. The neutron stars 
will no longer be causally disconnected from our dynam- 
ical boundary conditions as was the case in the short 
timescale simulations performed in the previous section. 
The computational boundary will therefore have a much 
greater affect on the simulation results if put at the same 
spatial location. One may need to greatly increase the 
distance from the center of mass of the system to the 
computational boundary, but one must be careful not to 
sacrifice the spatial resolution at the same time. 

Preliminary results obtained in [s^ suggest that wave- 
forms calculated from a numerical relativity binary inspi- 
ral simulation can be done in an accurate fashion (e.g., 
with errors less than 10%) only when the extraction ra- 
dius is approximately one gravitational wavelength Xg 
away from the center of mass of the system (for compar- 
ison, the NS-4 calculations from the previous section has 
the outer boundary at only O.OSAgu,). However, gravi- 
tational wave extraction is just one aspect of our sim- 
ulation requirements. Other effects including spacetime 
dynamics and dynamics of the binary may require the 



Configuration 


grid size 


Ax /Mo 






NS-A 


643x643x325 


0.2085 


0.257 


0.257 


NS-B 


323x323x165 


0.2085 


0.129 


0.257 


NS-C 


313x313x160 


0.2607 


0.156 


0.160 


NS-D 


259x259x133 


0.2607 


0.129 


0.160 


NS-E 


163x163x85 


0.4171 


0.129 


0.160 



TABLE II: The computational domain configurations used 
for the large timescale binary neutron star general relativis- 
tic simulations. All large timescale simulations are performed 
with a CFQE configuration characterized by an orbital an- 
gular velocity of Q.Mo — 0.01204, where the geodesic sepa- 
ration of the neutron stars is ^i_2 = 27.57Mo. The gravita- 
tional wavelength Xg^ corresponding to this configuration is 
Xgw = 5(17) = 260.9Afo. Tft denotes the (coordinate) dis- 
tance from the center of the orbiting binary to the bound- 
ary of the computational domain, denotes the (coordi- 
nate) distance between the center of the orbiting binary to 
the boundary of the computational domain used in solving 
for the CFQE initial data configuration. 



boundary to be even farther away. These facts, coupled 
with the complexities involved in solving the full Einstein 
field equations by computer, render the problem of ob- 
taining simulations accurate enough to probe the details 
of large-separation orbiting binary neutron stars a most 
difficult one. 

In this section we analyze the numerical evolutions of 
one particular CFQE initial configuration which has a 
larger geodesic separation £1^2 than that of the CFQE 
ISCO configuration. Specifically, the angular orbital ve- 
locity of the CFQE configuration we use exclusively in 
this section is ilMg = 0.01204. This configuration has a 
geodesic separation of ^1,2 = 27.57Afo, and corresponds 
to the second smallest data point shown in Figure 
According to the study in section Hvl this configuration 
has a violation of the QE and CF assumptions at the 
22% and 3% levels, respectively, soon after the evolution 
starts. We numerically evolve this CFQE initial data 
configuration using our fully consistent general relativis- 
tic treatment. The gauge conditions used for these simu- 
lations are the "l-flog" equation for the lapse a fEa. l35() 
and Eq. 1201 for the shift vector /3* . There is no need to 
use maximal slicing in this section as comparing to the 
CFQE-sequence is no longer the point. 

In Table im we list the properties of the various compu- 
tational domains, varying both the resolution and outer 
boundary placement, used for our long timescale numer- 
ical evolutions. Our numerical implementation allows us 
to use a different location for the outer boundary of the 
computational domain for our initial data solve of the 
CFQE configuration as that used for the dynamical evo- 
lution. We denote r^^ as the shortest coordinate distance 
between the center of our computational domain and the 
computational boundary of our cubical domain used in 
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FIG. 17: We plot the evolution of the geodesic distance be- 
tween the maximum rest mass density of the two neutron 
stars as a function of coordinate time for various resolutions 
and outer boundary placements (see TablelTTlfor configuration 
specifications) . 



solving the initial data problem for the CFQE configu- 
ration. We denote rb as the shortest coordinate distance 
between the center of our computational domain and the 
computational boundary of our cubical domain used in 
our fully consistent general relativistic numerical simula- 
tions. Note that in Table ITU the number of grid points 
refers to that used in the full dynamical evolution. For 
those computational domains where r^d > r?,, we have 
used the same resolution Aa; with a larger number of grid 
points to solve for the CFQE initial data configuration. 

In Figure El we plot the geodesic separation £1^2 as 
a function of coordinate time t using this initial data 
for the various computational domains listed in Table ITTI 
Qualitatively, we observe a number of interesting fea- 
tures. First of all, note that all of the simulations in 
Figure ITTI displav an eccentricity in the orbit. Some level 
of eccentricity is expected due to the fact that the CFQE 
configuration is constructed with the explicit assumption 
that the time derivative of the separation exactly van- 
ishes, as argued in [13, 13- That is, the assumption of 
the radial velocity of the binaries to be exactly zero in- 
stantaneously (which is done by assuming the existence 
of a timelike helical Killing vector field in the construc- 
tion of CFQE configurations) is not consistent with the 
astrophysically relevant scenario of quasicircular binary 
evolution, where the magnitude of the radial velocity, 
while small, never vanishes. One can obtain some idea re- 
garding the level of orbital eccentricity intrinsic to CFQE 
configurations by analyzing the point-particle dynamics 
in the post-Newtonian approximation (54] . In |54l |. it is 
shown in the context of the post^/^-Newtonian approxi- 
mation that if the assumption of circular motion is used 
to construct the initial conditions (that is, the conditions 
f = and f = are used to specify the initial conditions) , 
then the resulting orbit will have a nonzero eccentricity 



whose value depends on the initial separation. 

Ideally, we would like to directly compare our general 
relativistic results for the eccentricity of the orbits appar- 
ent in Figure ITTI with the post-Newtonian results in [5J|. 
The reason we are unable to do this at the present time 
is obvious from Figure ITTI Note the large difference be- 
tween the eccentricity for computational domain config- 
uration NS-A as compared to the other configurations. 
In fact, it is obvious from Figure El that the NS-A con- 
figuration is qualitatively different from the other con- 
figurations. Recall that all configurations NS-A through 
NS-E use the same CFQE configuration as initial data; 
note from Table UTI that the only parameters that differ 
in producing the results from Figure 1171 are the spatial 
resolution of the discretized computational domain, the 
location of the outer boundary of the computational do- 
main during numerical evolution where the dynamical 
boundary conditions are imposed, and the location of 
the outer boundary of the computational domain used 
in solving for the CFQE initial data configuration. A 
cursory study of Table |n| reveals the major difference 
between computational domain configuration NS-A and 
the others, namely, the location of the outer boundary 
Tft of the computational domain during numerical evo- 
lution. The boundary distance of the computational 
domain configuration NS-A is slightly larger than 1/4 
of the gravitational wavelength Xg^ characterized by the 
CFQE initial data configuration, while all of the other 
configurations have rf, < 0.156 \gw 

We must always quantify the errors of any numerical 
calculation. Here, our numerical errors originate from 
two distinct sources. The first source of error is the trun- 
cation error associated with our discretization parame- 
ter Ax, while the second source of error is induced by 
the outer boundary conditions imposed on our numer- 
ical calculation. While the behavior of the truncation 
error is well understood (it is a local error which scales 
as an integer power of the discretization parameter Ax), 
the assessment of the effect of the boundary conditions 
is not as straightforward. Theoretically, one could place 
the outer boundary of the computational domain suffi- 
ciently far away so that the outer boundary condition 
would not be causally connected to the compact objects, 
and thus would have no effect on them during numeri- 
cal evolution. However, this is not a practical solution, 
due to the limitations of computational resources (espe- 
cially for a unigrid code; adaptive mesh refinement could 
be used in this direction). As our outer boundary is 
causally connected to the neutron stars in our simula- 
tions, we must attempt to assess the errors introduced 
by the outer boundary conditions in our numerical sim- 
ulations. We assume that the error induced by the outer 
boundary conditions can be expanded in terms of pow- 
ers of l/rb, and that the error goes to as r;, — > 0. We 
therefore assume an error function for the eccentricity e 
as 

e„ = ee.act + Ci(Ax) -I- C2(Ax)' -f ^ + (77) 



where e„ denotes the measured value of eccentricity from 
our numerical solution using discretization parameter 
Ax and outer boundary location r^. Using the defini- 
tion of eccentricity defined in '53| in which the eccen- 
tricity of the orbit is calculated from the orbital sep- 
aration as a function of time, we compute the eccen- 
tricity associated with each simulation shown in Fig- 
ure E] We find that {en)j^s-A = 0.0124, (e„)Ars_B = 



0.0327, (e„) 
(e«) 



NS-C 



NS-E 



0.0397, (e„)^s_c = 0.0434, and 
0.0605. We can then solve for the unknown 
quantities Cexact, Ci, C2, C3, and C4 in Eq. [23 We 
find that the Richardson extrapolated value of the ec- 
centricity is Cexact = —0.127, and that the leading error 
terms for the truncation error and boundary error are 
CiAx — 0.067 and Cs/ri, = 0.11, respectively. Above 
and beyond the fact that the Richardson extrapolated 
value of the eccentricity Cexact is negative, the obvious 
sign that we are not in the convergence regime (i.e. that 
the higher order terms neglected in the error expansion 
of Eq. [77| are not relatively small) is that \eexact — Gn\ 
is larger than the error terms in Eq. 23 As the com- 
putational resources available at the present time do not 
currently allow us to use our unigrid code to simultane- 
ously decrease the discretization parameter Ax further 
and increase the distance rt, from the center of mass to 
the outer boundary, we must admit that we can, at this 
time, make no definite conclusion as to the inherent ec- 
centricity in CFQE configurations used as initial data in 
numerical relativity. However, the prospect of being able 
to determine this point in the near future is good; we 
can expect both mesh refinement techniques and better 
outer boundary conditions to greatly aid in reducing the 
errors induced by the outer boundary in our numerical 
calculations. 



A. orbital decay rate 



Recall from Sections IIII Bl and IIII CI that the binding 
energies Ei, and fEas.l58land l63l respectively) shown 
in Figures and El represent the binding energy of each 
neutron star binary CFQE configuration as a function 
of geodesic separation ^1^2 (the orbital angular velocity 
r2 is monotonically increasing with decreasing geodesic 
separation ^1,2)- In the CFQE-sequence approximation, 
this binding energy is slowly converted to gravitational 
wave energy, and Figures [S] and tell us how much 
gravitational radiation energy is produced for changes 
in geodesic separation. We can approximate the rate of 
energy loss at any specific point in this sequence using 
the standard quadrupole formula (see, e.g., |39|). which 
reduces to 
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(78) 



for two point particles of mass M in circular orbit with ra- 
dius R and orbital angular velocity Vl. We interpolate the 
data represented in Figures [3 and |^ with a cubic spline 
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FIG. 18: The geodesic separation £1,2 as a function of time. 
Shown are the results from our fully consistent general rela- 
tivistic calculation NS-A and NS-B (see Table^J ■ Also shown 
are predictions from the CFQE-sequence approximation, i.e. 
the curve obtained by numerically integrating Eq. 1791 The 
curve labeled "CFQE sequence" was constructed using the 
standard definition of binding energy, Eq. 1581 (see Figure |3 , 
and is terminated at the CFQE ISCO point (£i,2/Mo = 24.0). 
The curve labeled "CFQE sequence [ASns ~ 0)" was con- 
structed using the binding energy _Ej, Eq. 1631 where the 
neutron star spin remains constant throughout the entire 
CFQE sequence (see Figure El a^ud the discussion in Sec- 
tion I^^Cl) ^nd is terminated at the neutron star touching 
point {li,2/Mo = 22.6). 



to obtain the effective binding energies Ef, (Eq. I58|) and 
E'l, (Ea. l63|l as a function of geodesic separation ii^2- We 
can then easily find dEi,/d£i,2 as a function of geodesic 
separation. An estimate of the time rate of change of 
geodesic separations is then 



d£i^2 dEg^/dt 



dt 



dEu 



'1,2 



(79) 



which can be numerically integrated to produce the 
geodesic separation ^1^2 as a function of time predicted 
from the CFQE-sequence approximation. We plot these 
functions in Figure ITHl along with results from our fully 
consistent general relativistic calculation NS-A and NS- 
B. 

Note that the standard CFQE-sequence prediction of 
the evolution of the geodesic separation £1,2, labeled 
"CFQE sequence" in Figure terminates at a neu- 
tron star geodesic separation £1^2 = 24.0 Mq, whereas 
the modified CFQE-sequence prediction, labeled "CFQE 
sequence [ASns = 0)", terminates at a smaller geodesic 
separation of €1.2 = 22.6 A^q. This is due to the fact 
that the standard effective binding energy Ei, defined by 
Eq. 1581 actually attains a minimum before the neutron 
stars touch (see Figure ISJ , whereas the effective bind- 
ing energy for the CFQE-sequence where the spin of the 
neutron stars do not change {E[^ defined by Eq. I63|l is 
monotonically decreasing as the neutron star separation 
decreases, right through to the point where the neutron 
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stars are touching (see Figure 13 • Therefore, the first 
"CFQE-sequence" prediction in Figure ITHl terminates at 
this minimum point (defined as the ISCO configuration 
of the sequence), whereas the second "CFQE sequence 
{ASns = 0)" prediction in the figure terminates at the 
CFQE configuration where the neutron stars are touch- 
ing. One might expect that a full solution to the Ein- 
stein equations using a CFQE configuration as initial 
data might, in fact, have an evolution of the geodesic 
separation of the neutron star binaries £1^2 that would lie 
somewhere in between the two CFQE-sequence approxi- 
mations shown in Figure [TSI After all, these two CFQE- 
sequence approximations represent the extreme cases of 
the evolution of the neutron star spin; the first represents 
complete tidal locking during the entire evolution, the 
second represents absolutely no change in the spin state 
of the individual neutron stars during the entire evolu- 
tion. Of course, the actual solution to the full Einstein 
equations would be expected to couple, to some amount, 
the increasing orbital angular frequency of the binary 
neutron stars to the spin of the individual stars. The 
magnitude of this coupling would determine which of the 
two CFQE-sequence approximations shown in Figure ITHl 
would be considered to be "more" correct. It is exactly 
this question which could be answered with our fully con- 
sistent general relativistic calculations. However, in Fig- 
ure ^1 we once again see the dramatic effect of the outer 
boundary conditions on our numerical simulations. Note 
that computational domain configurations used in the 
calculations NS-A and NS-B have the same spatial res- 
olution Ax = 0.2085 Mq and outer boundary location 
rid = 0.257 Xgw used for solving the CFQE initial data 
configuration. The only difference between the configu- 
rations NS-A and NS-B is the outer boundary placement 
used during the dynamical evolution: the outer boundary 
of the computational domain used during numerical sim- 
ulation NS-A is twice as far away from the center of mass 
as that used in simulation NS-B. One may be tempted to 
conclude that the calculation NS-A could be said to be 
"more correct" than NS-B, since the outer boundary is 
located farther away for the NS-A configuration. How- 
ever, one must be careful when trying to apply physical 
intuition to numerical results that are not in the con- 
vergence regime. In this case, for instance, there is no 
reason to believe that the error induced by the outer 
boundary on our calculation is a monotonic function of 
rt (whereas, e.g., it is true in general that the trunca- 
tion error is monotonic in Ax as Ax — *■ 0). Due to the 
wave nature of the gravitational radiation being emitted 
by the binary neutron stars, the boundary errors could 
have an oscillating component. This makes it even more 
difhcult to try to do a Richardson extrapolation type of 
error analysis in realistic compact object simulations in 
numerical relativity. 

While it should be expected that the extraction of 
gravitational radiation for a numerically generated space- 
time could be highly sensitive to the location of the outer 
boundary (see 38]), there have been, up to now, no re- 



sults showing what effect the outer boundary conditions 
can have on the details of the orbits of compact binaries 
on timescales larger than one orbital period. Here, we see 
that the dynamical outer boundary conditions can, and 
do, significantly affect the orbital parameters of compact 
binaries during numerical evolution. 



VI. CONCLUSIONS 

To date, the only fully general relativistic simulations 
of corotating binary, F = 2 polytropes is js^l • The study 
in 37] differs from the present study in several impor- 
tant ways. One basic difference is that while we focus on 
the capability of simulating astrophysical realistic neu- 
tron star binaries, jSTl ] focuses on the dynamics of the 
final merger of the two neutron stars. Thus, the initial 
data used in [s^ are CFQE configurations either at ISCO 
separation, or closer than ISCO separation (ISCO sepa- 
ration here means simply the neutron star separation of 
the unique configuration that corresponds to a minimum 
of the binding energy in the constant rest mass CFQE- 
sequence, e.g., in Figure|3J|. Also, in the study [s^], it was 
found necessary to manually decrease the orbital angular 
momentum of each CFQE configuration at the level of 
several percent, in order to precipitate the binary merger 
more quickly. Complementary to that approach, we are 
trying to use numerical relativity as a tool to assess the 
fidelity of the CFQE-sequence approximation itself. We 
found that in order to perform simulations of the neutron 
star binary systems compatible with realistic astrophysi- 
cal scenarios, one must perform simulations using initial 
data at a distance considerably larger than the ISCO 
separation when using corotating CFQE configurations 
as initial data. 

We have outlined a generic method for analyzing the 
regime of validity of the CFQE-sequence approximation 
and have applied this method to the case of equal mass, 
corotating binary neutron stars. We have found that, 
for corotating neutron stars, the violation of the time- 
like helical Killing vector field existence assumption was 
an order of magnitude larger than the violation of the 
assumption of conformal flatness. Specifically, we have 
demonstrated that initial data specified by a CFQE con- 
figuration with neutron stars having an initial geodesic 
separation of less than 47 Mq (which is slightly more than 
6 neutron star radii, or roughly twice the ISCO configu- 
ration separation) would produce a solution to the Ein- 
stein field equations that violates the Killing vector field 
assumption by more than I part in 10; the conformal 
flatness assumption would be violated by more than I 
part in 100. We thus conclude that, in the corotating 
case, the CFQE-sequence approximation for neutron star 
separations of 47 Mq and less violates the Einstein field 
equations at a level larger than 10%, and thus numerical 
simulations starting with similar CFQE configurations as 
initial data cannot, therefore, be considered as approxi- 
mating a realistic neutron star binary inspiral. 
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We note that the violations of the assumptions of the 
CFQE-sequence approximation that we observe in our 
general relativistic calculations for the corotating binary 
systems occur on timescales that are two orders of mag- 
nitude shorter than the orbital timescale. We suspect 
that this may be due to interactions between the spin as- 
sumption of the individual neutron stars (the corotating 
assumption) and the CFQE assumptions. We have shown 
that the characteristic shape of the effective binding en- 
ergy curve within the CFQE approximation is highly sen- 
sitive to the spin kinetic energy of the individual neutron 
stars. We have shown that if we subtract out the spin 
kinetic energy of the neutron stars in the construction 
of the effective binding energy (which approximates the 
case where the neutron star spin does not increase as 
the orbital angular velocity increases), then the resulting 
binding energy curve will have no minimum, and thus 
the CFQE-sequence approximation would not predict the 
existence of an ISCO configuration. We speculate that 
specifying neutron stars with irrotational spin states in 
the CFQE-sequence approximation may yield a smaller 
violation of the Einstein field equations for a fixed neu- 
tron star separation. The analysis we have developed in 
this paper can be used for a detailed investigation of this 
effect. More interestingly, the analysis we have developed 
might provide a way to determine a spin state most con- 
sistent with the CFQE approximations, and hence pro- 
vide a more realistic set of initial data that can be used 
to start simulations at a smaller initial separations. 

We have shown that, for our specific neutron star mod- 
els, we require a resolution of approximately Ax — O.IA/q 
in order to adequately resolve the neutron stars ("ade- 
quately resolve" here refers to verifying that we are in the 
convergence regime through an appropriate Richardson 
extrapolation technique; this is a much more stringent 
condition than has been typically used in numerical rel- 
ativity studies to date involving neutron stars and black 
holes) . This resolution scale is over three orders of mag- 
nitude smaller than the characteristic wavelength of the 
gravitational radiation emitted during the last five to ten 
orbits of the neutron star inspiral process. 

We have also shown that the location of the outer 
boundary of the computational domain can have a signif- 
icant impact on the details of the evolution of the com- 
pact objects on timescales of the orbital period. Specif- 
ically, we have seen that changing the linear dimensions 
of our computational domain from 0.3 Xgw to 0.5 Xg^^ 
can significantly impact the dynamics of binary neutron 
stars during the first several orbits. This should serve as 
a warning to the numerical relativity community study- 
ing simulations of compact binaries with the hope of ex- 
tracting gravitational wave information: not only will the 
outer boundary inhibit the actual process of extracting 
the gravitational waves, but they also directly affect the 
sources of the gravitational waves themselves. While our 



dynamical boundary conditions are not the best choice, 
and a more consistent treatment, e.g. constraint pre- 
serving boundary conditions, would most likely improve 
the situation, it may be that numerical relativists will 
be forced to push the outer boundary of the computa- 
tional domain to the "local wave zone" (which in this 
case means rf, > Xgw) in order to provide realistic gravi- 
tational waveforms suitable for use as templates in grav- 
itational wave detectors. 

Unfortunately, this makes the numerical study of or- 
biting compact objects in numerical relativity particu- 
larly hard. Note that every time we increase the reso- 
lution in our 3D simulations by a factor of two, keeping 
the outer boundary fixed, we must use 16 times the 
amount of computational resources to perform any par- 
ticular simulation. Also, every time we increase the outer 
boundary distance by a factor of two, we must use 8 
times the amount of computational resources. Therefore, 
if we wanted to simultaneously increase both the resolu- 
tion and the outer boundary distance by a factor of 2, we 
would require over two orders of magnitude more com- 
putational resources. While wc have shown in this paper 
that it is possible to track the details of finite sized com- 
pact objects in full numerical relativity, what remains is 
to be able to do so in such a way that all of the numeri- 
cal errors (both truncation and boundary errors) can be 
demonstrated to be small over a timescale of several or- 
bital periods. This will be an extremely challenging task, 
given the current level of computational resources avail- 
able. It may be necessary to employ mesh refinement in 
order to accurately simulate all of the physical degrees of 
freedom that we are interested in. 
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